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We present a calculation of the charmonium spectrum with three flavors of dynamical staggered quarks from 
gauge configurations that were generated by the MILC collaboration. We use the Fermilab action for the valence 
charm quarks. Our calculation of the spin-averaged IP-IS and 2S-1S splittings yields a determination of the 
strong coupling, with a^g(Mz) = 0.119(4). 



1. INTRODUCTION 

The current experimental program of precision 
flavor physics at the B factories, at CESR-c, and 
at the Tevatron needs accurate lattice QCD cal- 
culations of the relevant hadronic matrix elements 
to yield stringent constraints on the CKM sector 
of the standard model. Precision lattice QCD re- 
sults in turn require that the systematic errors as- 
sociated with lattice calculations be brought un- 
der control to the desired accuracy. The most im- 
portant sources of systematic error in lattice cal- 
culations include the incomplete inclusion of sea 
quarks (quenched approximation) , lattice spacing 
artifacts, perturbative errors, and the chiral ex- 
trapolation. We are planning a series of lattice 
calculations of the phenomenologically most im- 
portant quantities in the B, D, cc and bb sys- 
tems. We will address the first two sources of 
systematic error by performing simulations with 
highly improved actions on gauge configurations 
with rif = 2 + 1 improved staggered quarks [Q. 
This effort must be complemented by the cor- 
responding perturbative matching calculations of 
the improvement coefficients and currrent renor- 
malizations. This should be possible with recent 
advances in automated perturbation theory . 



*Talk and poster presented by P. Mackenzie and D. Men- 
scher. 



The MILC collaboration has generated dynam- 
ical gauge configurations Q using improved stag- 
gered and gluon actions. Their configurations in- 
clude three (or 2 + 1) flavors of light staggered 
fermions at several different light quark masses 
ranging from m s to m s /5. Hence, the system- 
atic errors usually associated with the quenched 
approximation should be absent with these con- 
figurations. Furthermore, since numerical simu- 
lations with rather light quark masses are feasible 
with staggered actions, the issue of chiral extrap- 
olations may be carefully studied. 

The heavy quark action used in this work is 
based on Ref. (§. It is related to NRQCD, but 
uses the four component fields and operators of 
the Wilson action rather than the two compo- 
nent fields and operators of NRQCD. Similar to 
NRQCD, the space-like and time- like components 
of the operators are uncoupled, and the coeffi- 
cients of the operators are mass dependent. This 
action smoothly interpolates between an ordi- 
nary light quark action as am — > and NRQCD 
when am > 1, but is applicable at all values of 
am. Our formalism can be regarded as a sum- 
mation of terms of the form (am) p to all orders 
in the normalizations of operators, which is use- 
ful when m 3> Aqcd- To 0(a), our action uses 
the same operators as the clover action Q . Start- 
ing at 0(a 2 ), the operators are somewhat differ- 
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ent from those in the analogous light quark ac- 
tion. The reason is that a two-hop correction 
to the Wilson time derivative operator cannot 
be used because it introduces ghost states for 
heavy quarks. Its effects must instead be dupli- 
cated with Hamiltonian-style operators ||. The 
action may be particularly useful for the charm 
quark on lattice spacings with am c < 1. It has 
smaller discretization errors than standard light 
quark actions since the (am c ) p errors, which it 
resums and eliminates, are much larger than the 
remaining (aAqcrj) p discretization errors. Since 
it has a well-defined a — > limit, it may be more 
convergent for am < 1 than NRQCD, which does 
not. 

Now that we are calculating with dynamical 
quarks, lattice spacings obtained from the sim- 
plest quantities should all agree within errors. We 
start our work with a study of the charmonium 
spectrum, to test our methods and lattices. As a 
byproduct, we obtain a new determination of the 
strong coupling. 

2. DETAILS OF THE CALCULATION 

MILC uses the Asqtad action [|j for the stag- 
gered fermions which contains errors of 0(a s a 2 ), 
and an improved gluon action with 0(p? s a 2 ) 
errors. For the charm quarks we use the 
Sheikholeslami-Wohlert action Q (which has 
0(a s a) errors) with the Fermilab interpretation 
Q. The quenched gauge configurations are gen- 
erated with the Wilson gauge action. Table |l| 
lists the simulation parameters for all three lat- 
tices. The quenched lattice has a slightly smaller 
lattice spacing than the MILC lattices. 

As usual, we calculate charmonium two-point 
functions using smeared source and sink op- 
erators. For this purpose, when working on 
MILC lattices, we use the Richardson poten- 
tial model wave functions shown in Figure |l|. 
The quenched propagators were generated using 
Coulomb wave functions. 

2.1. Fits 

As shown in Refs. J8|||, constrained fits allow 
for better control of the systematic error due to 
excited state contributions, because they allow us 
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Simulation parameters for the three lattices. 



Size 


16 3 x 32 


20 3 x 64 


20 a x 64 


71/ 





3 


2+1 





5.9 


6.85 


6.76 


configs 


300 


174 


298 


am s 


00 


0.05 


0.05 


ami 


00 


0.05 


0.01 




0.1227 


0.113 


0.113 


wf's 


5, IS, 2S 


S, IS 


5, IS 



to fit the correlators to a large number of states 
without loss of accuracy. Furthermore, with con- 
strained fits one is also able to use all of the time 
slices in the fits without adjusting t m m and t max 
to determine the best fit. We investigate this is- 
sue by comparing constrained and unconstrained 
fits. We fit meson correlators to the form 

G(t-{Z n },{E n })= (1) 
£^(e-^ + e -^-')). 

n 

We force our energy levels to be ordered by defin- 
ing 

AE n = E n ~ £„_i = exp(e„) . (2) 
In our fits we use Baysian statistics, 

2 2 _ 2 , 2 lo\ 

X * Xaug — X ~r Xprior ! \*'J 

where Xprior i s use d to constrain Eq, e„, and Z n 
to a predetermined range. 

For simplicity, we test our fit program on the 
quenched lattice. Figures | | show our fit results 
for the rj c ground and first excited state energies 
as functions of the number of states in the fit (n) 
obtained from the r] c propagator with a (5-function 
source and sink. We observe that the fit results 
- particularly the error bars — for the ground 
and first excited state energies are stable under 
adding more states to the fit when constrained 
fits are used, but not in the case of unconstrained 
fits. The unconstrained fits do not use Eq. ||, 
but they do use the energy ordering constraint 
of Eq. ||, which is probably the reason for the 
relative stability of the central fit values even in 
the unconstrained case. 
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Figure 1. Wavefunctions for the 15*, 2S, and IP 
states in charmonium. 



The prior constraints must be chosen so that 
they do not unduly influence the physical results. 
Figure || shows the dependence of the ground 
state fit on the prior width. The fit includes four 
states and the ranges for all priors are varied to- 
gether with the ground state energy prior width. 
We see that the fit results are stable, once the 
prior width is large enough. Furthermore, the er- 
ror on the ground state energy is unaffected by 
the variation of the prior width, after the plateau 
is reached. Table |] lists typical choices for the 
energy prior values and ranges used in our fits. 
The prior values for the Z n are chosen by match- 



Table 2 



Typical 


energy prior 


values and 


ranges. 


n f 





3 


2 + 1 


rjc- E 


1.4 ±0.2 


1.9 ±0.2 


1.9 ±0.2 


AE„ 


0.4 ±0.2 


0.4 ±0.2 


0.4 ±0.2 


h c - Eq 


1.7 ±0.2 


2.3 ±0.2 


2.3 ±0.2 


AE n 


0.4 ±0.2 


0.5 ±0.2 


0.4 ±0.2 



Figure 2. Comparison of fit results from con- 
strained and unconstrained fits. Shown is the 
rj c ground state energy — obtained from the local- 
local correlator as a function of the number of 
states in the fit (n). o: constrained fit, x: uncon- 
strained fit 



ing the hadron propagators evaluated at t = 1 to 
Eq. [I]; the range is usually set to a factor of three 
of the central value. 

The results discussed in the next section are 
obtained from fits to multiple correlators, mak- 
ing use of the different source and sink operators 
listed in Table |l|. These fits are generally con- 
sistent with fits to the delta-function correlators, 
albeit with smaller statistical errors. 

3. THE SPECTRUM 

Figure ^| summarizes our results for the char- 
monium spectrum. Our results for the hyperfine, 
IP-IS, and 2S-1S splittings are shown in Fig- 
ures as functions of the light quark mass. We 
observe very little light quark mass dependence 
in the hyperfine splitting, and our chirally ex- 
trapolated result still disagrees with experiment, 
although the inclusion of dynamical quarks has 
removed one possible cause of a small hyperfine 
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Figure 3. Same as Figure g, but for the first ex- 
cited state, the rf c . 



Figure 4. Variation of the r\ c ground state energy 
fit result with the prior width, i is defined via 
a = 2Vo, where a and <tq are the varied and 
standard choice of prior widths respectively. 



splitting: the fact that the short distance coupling 
constant is too small in the quenched approxima- 
tion. We note that the charm quark action is only 
0(a) improved. The leading order operator which 
controls the spin splitting is ^po^vF^ip. Its coef- 
ficient is being included only at tadpole improved 
tree-level, and it is plausible that the observed 
discrepancy is a result of both 0{a s a) errors and 
0(a 2 ) lattice spacing artifacts. We will be able to 
study this issue further once our 0(a 2 ) improved 
action Q is ready for numerical simulations. 

The spin-averaged IP-IS and 2S-1S splittings 
are considerably less sensitive to the leading or- 
der lattice artifacts; they are used to determine 
the lattice spacing. Figures 0-^ indicate that the 
dependence of the spin-averaged splittings on the 
light quark mass is mild. 

The extent to which the lattice spacings from 
the IP-IS and 2S-1S splittings disagree with each 
other is an indication of residual systematic er- 
rors in our simulation, which are a combination 
of higher order lattice spacing and sea quark ef- 
fects. For this comparison we compute the ratio 



of lattice spacings, 



lal 



R 



AM(2S-1S) 
_ AM(IP-IS) 
= AA/(2S-1S) cxp 

AM(IP-IS) 



(4) 



Table || compares the lattice spacings obtained 
on all three lattices. On the quenched lattice the 
deviation of R from unity is about 1.5 standard 
deviations. On the two MILC lattices, the devi- 
ation of R from unity is less significant, because 
of the still somewhat large statistical errors. To 



Table 3 

The lattice spacings from the spin-averaged IP- 
IS and 2S-1S splittings for the three lattices with 
statistical error bars. 



n f 


a-^lP-lS) 
(GeV) 


a- 1 (2S-lS) 
(GeV) 


R 




3 

2 + 1 


1 7OQ+0.062 

1. 1 O0_Q Q3Q 

1 70+ 016 

1 -' u -0.18 
l.OD^_ 053 


1 RQ7+U.UV9 
1 - UJI -0.059 
1 co+0.15 
1 - oz -0.11 

i 42fi+ 073 


-, i 9 +0.14 
i - LZ -0.17 

i 0Q7+ - 073 

J-.USM -0.068 
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Figure 5. The charmonium spectrum in compar- 
ison, o: rif = 0, o: = 3, □: n./ = 2 + 1. 



140 
120 
| 100 

a so 

60 



0.00 



0.05 
am, 



0.10 



Figure 6. The hyperfine splitting vs. ami (cir- 
cles)). Shown also is the quenched result (square), 
positioned in the plot at finite ami for illustration, 
and the experimental result (burst) positioned at 
ami = 0. 
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Figure 7. The IP-IS splitting vs. ami (circles). 
Shown also is the quenched result (square), posi- 
tioned in the plot at finite ami for illustration. 
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Figure 8. The 2S-1S splitting vs. ami (circles). 
Shown also is the quenched result (square), posi- 
tioned in the plot at finite ami for illustration. 



clarify the situation we need to reduce the statis- 
tical errors of the results on the MILC lattices. 

Now that dynamical fermions are included in 
the calculations, lattice spacings from the best- 
understood quantities should be consistent. We 
note that although the lattice spacings obtained 
from the IP-IS and the 2S-1S splittings are con- 
sistent with each other, the lattice spacing from 
the IP-IS splitting is consistent with several de- 
terminations of the lattice spacing from the T 
system jlQ] , while the one obtained from the 2S- 
1S splitting is not. It is possible that higher or- 
der discretization effects are responsible for this. 



A more interesting possibility is that due to the 
fact that the 2S is so close to the DD threshold, 
physical effects from the coupling to DD channels 
may have a more dramatic effect in the 2S than 
in other states. 

4. THE STRONG COUPLING 

The spin-averaged splittings discussed in the 
previous section are used to determine the strong 
coupling a s . The 2 + 1 flavor lattices are our 
most realistic. We take the IP- IS splitting as our 
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most reliable determination of the lattice spacing 
because of the possible threshold effects in the 2S 
state. Following the procedure of Ref. we 
obtain the strong coupling from the plaquette. 
For our actions we have Q 

- ln(Tr Up) = (5) 
3.0682 a P {q*) [1 - a P (0.770 + 0.09681 n f )] 

where q* = 3.33/a is the BLM scale for nr = 3. 
The coupling ap is defined to coincide through 
one-loop order with ay, the coupling defined from 
the heavy quark potential. Using 



2 

1 H — ap 

TT 



0{a 2 P ) 



we obtain 



a m (M z ) =0.119 ±0.004. 



(6) 



(7) 



The difference between our value and the value re- 
ported in Ref. [ O] arises mainly from differing im- 
plicit treatment of 0(a 3 ) corrections, which will 
soon be known. The main sources of uncertainty 
are 0(a 3 ) corrections (3%), discretization errors 
(2%), and statistical errors (1%). All three should 
be significantly reduced soon. 

5. CONCLUSIONS AND OUTLOOK 

We present preliminary results of a calculation 
of the charmonium spectrum on gauge configura- 
tions generated by the MILC collaboration using 
0(a 2 ) improved actions for the gluons and the 
nf = 2 + 1 dynamical staggered fermions. We use 
the O(a) improved clover action with the Fermi- 
lab interpretation for the charm valence quarks. 
Since the MILC configurations were generated 
with the correct number of sea quarks, no ex- 
trapolation in nf is necessary. Furthermore, the 
nondegenerate strange and light quark masses in 
the MILC lattices allow us to consider the chiral 
limit. Comparing results at am s = ami = 0.05 
and am s = 0.05, ami = 0.01, we find only mild 
light quark mass dependence for all spectral quan- 
tities we consider. Finally, our calculation yields 
a new determination of the strong coupling where 
systematic errors due to sea quark effects are un- 
der control. 



For future work, we are planning to improve the 
statistical accuracy of this work. The improve- 
ment of the heavy quark action beyond O(a) is in 
progress ||. 
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